/***************************************************************************
 *
 *   Copyright (C) 2005 - 2017 by Willem van Straten
 *   Licensed under the Academic Free License version 2.1
 *
 ***************************************************************************/

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include "Pulsar/DeltaRM.h"
#include "Pulsar/DeltaPA.h"

#include "Pulsar/Archive.h"
#include "Pulsar/Integration.h"
#include "Pulsar/PolnProfile.h"

#include <iostream>
using namespace std;

// #define _DEBUG 1

Pulsar::DeltaRM::DeltaRM ()
{
  threshold = 3.0;
  used_bins = 0;

  policy = new FrequencyIntegrate::EvenlyWeighted;
}

Pulsar::DeltaRM::~DeltaRM ()
{
}

//! Set the archive from which to derive the refined rotation measure
void Pulsar::DeltaRM::set_data (const Archive* archive)
{
  data = archive;
}

#if FIX_PGPLOT

#include "Pulsar/StokesCylindrical.h"
#include <cpgplot.h>

static void plot (const Pulsar::Archive* archive, isubint, unsigned chan)
{
  string name = "band" + tostring(chan) + ".ps/cps";

  cpgbeg (0, name.c_str(), 0, 0);
  cpgslw(2);
  cpgsvp (.1,.9, .1,.9);
  
  Pulsar::StokesCylindrical plot;

  plot.set_subint (subint);
  plot.set_chan (chan);
  plot.plot (archive);

  cpgend();
}

#endif

//! Set the phase bins to be included in the mean
void Pulsar::DeltaRM::set_include (const std::vector<unsigned>& bins)
{
  include_bins = bins;
}

//! Set the phase bins to be excluded exclude the mean
void Pulsar::DeltaRM::set_exclude (const std::vector<unsigned>& bins)
{
  exclude_bins = bins;
}

//! Set the policy used to integrate frequency channels
void Pulsar::DeltaRM::set_policy (FrequencyIntegrate::Divided* _policy)
{
  policy = _policy;
}

//! Refine the rotation measure estimate
void Pulsar::DeltaRM::refine ()
{
  if (!data)
    throw Error (InvalidState, "Pulsar::DeltaRM::refine", "no data");

  Reference::To<Integration> clone = data->get_Integration(isubint)->clone();

  FrequencyIntegrate fscr;
  
  policy->set_ndivide (2);
  fscr.set_range_policy (policy);
  fscr.transform (clone);

#if FIX_PGPLOT

  /* Following the creation of libpsrplot, plotting functionality will
     likely have to be added as a plugin; e.g. in the main() of the
     executable that uses this class */

  cerr << "Pulsar::DeltaRM::refine plotting half-band integrated profiles" << endl;

  plot (data, isubint, 0);
  plot (data, isubint, 1);

#endif

  /*
    remove the remaining Faraday rotation, so that residual delta_PA
    to be measured is (hopefully) less than 2pi 
  */
  clone->defaraday ();

  if (clone->get_nchan() != 2)
    throw Error (InvalidState, "Pulsar::DeltaRM::refine",
		 "nchan != 2 after FrequencyIntegrate");

  double f0 = clone->get_Profile(0,0)->get_centre_frequency();
  double f1 = clone->get_Profile(0,1)->get_centre_frequency();

  // speed of light in m/s
  double speed_of_light = 299792458;

  double lambda_0 = speed_of_light / (f0 * 1e6);
  double lambda_1 = speed_of_light / (f1 * 1e6);

  cerr << "frequency_0 = " << f0 << " MHz"
    "  ->  lambda_0 = " << lambda_0 << " m" << endl
       << "frequency_1 = " << f1 << " MHz"
    "  ->  lambda_1 = " << lambda_1 << " m" << endl;

  Reference::To<PolnProfile> profile0 = clone->new_PolnProfile(0);
  Reference::To<PolnProfile> profile1 = clone->new_PolnProfile(1);

  Pulsar::DeltaPA delta_pa;

  delta_pa.set_threshold (threshold);
  delta_pa.set_include (include_bins);
  delta_pa.set_exclude (exclude_bins);
  if (onpulse_weights)
    delta_pa.set_onpulse (onpulse_weights);
  
  Estimate<double> delta_PA = delta_pa.get (profile0, profile1);  
  Estimate<double> delta_RM = delta_PA/(lambda_1*lambda_1-lambda_0*lambda_0);

  rotation_measure = delta_RM + data->get_rotation_measure();
  
  cerr << "delta PA = <PA_1 - PA_0> = " << 180.0/M_PI * delta_PA << " deg.\n"
       << "delta RM = " << delta_RM << endl
       << "final RM = " << rotation_measure << endl;

  used_bins = delta_pa.get_used_bins();
}

